CRS Prep & Review

Setup

Data describing commercial landings (tonnes) in 2017.

# Read in commercial landings data
fish <- rast(here("Reference", "CRS", "commercial_landings_2017.tif"))

## Take a look
plot(fish)

# lets log so we can visualize this a little better
plot(log(fish+1))

How many fish were captured?

global(fish, "sum", na.rm = TRUE)

##                               sum
## commercial_landings_2017 86331119

# 86.3 million tonnes

Check against Sea Around Us to make sure this seems generally reasonable.

seems a little low https://www.seaaroundus.org/data/#/global?chart=catch-chart&dimension=sector&measure=tonnage&limit=10 but generally makes sense (close enough)
around 100 million https://www.seaaroundus.org/data/#/global?chart=catch-chart&dimension=sector&measure=tonnage&limit=10

What is the area of each cell?

cell_area <- cellSize(fish)
plot(cell_area, box = TRUE)

Describe what the plotted image is saying?

cell size is much smaller at the poles, increases as you head in towards the equator

Given this, what is kind of deceptive about these plots?

if density of fish landings in tonnes, pole values would appear higher than they are – apparent tonnes per cell would be higher
if count, same thing (seems higher than in reality)

plot(fish)

# lets log so we can visualize this a little better
plot(log(fish+1))

the cells in the northern and southern extremes have a much smaller area than those concentrated in the center
this may hide higher commercial fishing values in the center band, as those pixels will appear smaller when plotted
value associated w the cell is changing when you convert – resampling by diff method (nearest neighbor, linear interpolation, bilinear)
poles may appear to have higher values than they actually do – same value over smaller area = higher density value
more coarse resolution at the poles can lead to higher catch shown over a greater area than in reality
when we sum values across cells, the area of that cell is relevant when we do a global summation:
- From the documentation on cellSize()

Compute the area covered by individual raster cells. Computing the surface area of raster cells is especially relevant for longitude/latitude rasters.

But note that for both angular (longitude/latitude) and for planar (projected) coordinate reference systems raster cells sizes are generally not constant, unless you are using an equal-area coordinate reference system.

For planar CRSs, the area is therefore not computed based on the linear units of the coordinate reference system, but on the actual area by transforming cells to longitude/latitude. If you do not want that correction, you can use transform=FALSE or init(x, prod(res(x)))

Let’s convert these data to another coordinate reference system:

moll <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
fish_moll <- terra::project(fish, crs(moll))

# take a look
par(mar = c(5, 4, 4, 0)) # set margins: bottom, left, top, right
plot(fish_moll)

plot(log(fish_moll + 1), main = "Transformed CRS: Mollweide")

plot(log(fish+1), main = "Original CRS: WGS 84 EPSG 4326")

## Count the total tonnes
global(fish_moll, "sum", na.rm = TRUE)

##                                sum
## commercial_landings_2017 389271065

# 389.3 million

# calculate cell area
cell_area_moll <- cellSize(fish_moll)

#par(mai = c(1, 1, 3, 5)) # set margins: bottom, left, top, right
#par(oma = c(1, 1, 3, 4))
#par(mar = c(1,1,1,1))
plot(cell_area_moll, main = "Mollweide projection cell size")

plot(cell_area, main = "Original projection cell size")

389 million tonnes is crazy

(cell_area_moll$area)

## class       : SpatRaster 
## dimensions  : 764, 1546, 1  (nrow, ncol, nlyr)
## resolution  : 23330.19, 23330.19  (x, y)
## extent      : -18037447, 18031019, -8861637, 8962624  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        :      area 
## min value   :         0 
## max value   : 547291328

# min 0
# max is 547 million
cell_area$area

## class       : SpatRaster 
## dimensions  : 347, 720, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -180, 180, -85.5, 88  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## varname     : commercial_landings_2017 
## name        :       area 
## min value   :  122442386 
## max value   : 3077249667

# min: 122442386 (122 million)
# max: 3077249667 (3 billion)

Compare to the original tonnes? What is happening? How can you improve this? - original tons: 86331119 - Mollweide tons: 389271065 - the cell size area is now relatively consistent across the globe – old cell sizes have stretched etc., to fit a standardized cell size (I tried to extend the margins to see the full scalebar label values but struggled) - perform global calculations before reprojecting

could take raster of cell size, multiply by original (if density) raster, then sum values
^ divide if count
tonnes / area = density, then project
take new cell size, multiply by density raster
sum at end to check again – might need to reduce resolution, summary before projecting data

# get density (tonnes/area)
fish_density <- fish / cell_area

plot(log(fish_density))

# reproject to Mollweide
fish_density_moll <- terra::project(fish_density, crs(moll))
fish_density_moll

## class       : SpatRaster 
## dimensions  : 764, 1546, 1  (nrow, ncol, nlyr)
## resolution  : 23330.19, 23330.19  (x, y)
## extent      : -18037447, 18031019, -8861637, 8962624  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        : commercial_landings_2017 
## min value   :             0.000000e+00 
## max value   :             4.908811e-05

plot(log(fish_density_moll))

# get area of cells
new_cell_size <- cellSize(fish_density_moll)
plot(new_cell_size)

# multiply density by area to get count (tonnes)
new_tonnes <- fish_density_moll * new_cell_size
plot(log(new_tonnes + 1))

# sum to check
global(new_tonnes, "sum", na.rm = TRUE)

##                               sum
## commercial_landings_2017 85253920

# 85 million!

close enough :)
hopefully you don’t care about Fiji data too much (gets distorted with Mollweide projection)

#View(fish)

summary(fish)

##  commercial_landings_2017
##  Min.   :     0.00       
##  1st Qu.:     2.43       
##  Median :    22.25       
##  Mean   :   651.53       
##  3rd Qu.:    84.65       
##  Max.   :125112.30       
##  NA's   :47150

fish_polygons <- terra::as.polygons(fish)

fish_df <- terra::values(fish_polygons)

# Get the cell numbers
cells <- 1:terra::ncell(fish)

# Get the x and y coordinates for each cell
xy <- terra::xyFromCell(fish, cells)

# Extract the values
values <- terra::values(fish)

# Combine everything into a data frame
fish_df <- data.frame(
  x = xy[,1],
  y = xy[,2],
  value = values
)

fish_no_na <- fish_df %>% 
  na.omit()

See what happens when we make multiple CRS conversions:

# Robinson projection
rob <- "+proj=robin +datum=WGS84 +units=m +no_defs"

fish_moll_rob <- terra::project(fish_moll, crs(rob))
plot(log(fish_moll_rob + 1))

# now back to lat long:
fish_moll_rob_latlon <- project(fish_moll_rob, fish)

plot(log(fish_moll_rob_latlon + 1))

plot(log(fish + 1))

# check to see if it exactly matches the start
check <- fish_moll_rob_latlon - fish
plot(log(check + 1))

check

## class       : SpatRaster 
## dimensions  : 347, 720, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -180, 180, -85.5, 88  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## varname     : commercial_landings_2017 
## name        : commercial_landings_2017 
## min value   :                -66229.40 
## max value   :                 35487.52

What is this telling us?

the values have changed notably between projections
some are highly inflated (large positive numbers), other deflated (large negative in new projection)
adding error when we flip between resolutions
they’re all just approximations (by nature)

fish_moll_rob

## class       : SpatRaster 
## dimensions  : 772, 1542, 1  (nrow, ncol, nlyr)
## resolution  : 22047.33, 22047.33  (x, y)
## extent      : -17003404, 16993576, -8462847, 8557690  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        : commercial_landings_2017 
## min value   :                      0.0 
## max value   :                 132159.3

fish_moll_rob_latlon <- project(fish_moll_rob, fish) # match template (convert crs and resolution)
fish_moll_rob_latlon

## class       : SpatRaster 
## dimensions  : 347, 720, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -180, 180, -85.5, 88  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## varname     : commercial_landings_2017 
## name        : commercial_landings_2017 
## min value   :                      0.0 
## max value   :                 104592.5

fish_moll_rob_latlon_v2 <- project(fish_moll_rob, crs(fish)) # match crs (just convert crs)
fish_moll_rob_latlon_v2

## class       : SpatRaster 
## dimensions  : 750, 1553, 1  (nrow, ncol, nlyr)
## resolution  : 0.2316815, 0.2318308  (x, y)
## extent      : -180, 179.8014, -85.87307, 88  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : commercial_landings_2017 
## min value   :                      0.0 
## max value   :                 127646.3

CRS Prep & Review

Anna Ramji

2024